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Abstract. We describe a method for initializing characteristic evolutions of the Einstein 
equations using a linearized solution corresponding to purely outgoing radiation. This allows 
for a more consistent application of the characteristic (null cone) techniques for invariantly 
determining the gravitational radiation content of numerical simulations. In addition, we are 
able to identify the ingoing radiation contained in the characteristic initial data, as well as in 
the initial data of the 3+1 simulation. We find that each component leads to a small but long 
lasting (several hundred mass scales) transient in the measured outgoing gravitational waves. 
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1. Introduction 

It is well known in numerical relativity that current practice for the setting of initial 
data introduces spurious radiation into the system, in both the 3+1 and the characteristic 
approaches. The error in the initial data leads to an initial burst of spurious "junk" radiation 
that results from solving the constraint equations on a single hypersurface, without knowing 
the past history of the radiation content. Common practice regards the signal as physical only 
after it has settled down following the burst arising from the initial data solution. While it is 
straightforward to handle the junk radiation in this way, a more serious issue is whether the 
spurious radiation content of the initial data leads to longer-term transients in the wave signal. 
This question has been considered before, but previous work on the long-term effect of the 
initial data is limited E E] [3 S) . 

Characteristic extraction is a method of invariantly measuring the gravitational wave 
emission of an isolated source by transporting the data to null infinity (J'^) using a null 
formulation of the Einstein equations li51l6ll7l[8ll9l[T01. Initial data is needed on a null cone in 
the far field region, say r > lOOAf . Previous work has taken the simplistic approach of setting 
the null shear J = everywhere, although a recent investigation sets J by the condition that 
the Newman-Penrose quantity ipo — 110|. Setting shear-free initial data is not necessarily 
incorrect physically - for example in the Schwarschild geometry in natural coordinates J — 
everywhere, and it is possible to construct a radiating solution with J = everywhere at a 
specific time. However in the generic case, a radiating solution has J ^ 0, and imposing 
J = in effect means that the outgoing radiation implied by the boundary data must be 
matched at the initial time by ingoing radiation. In principle, the spurious incoming content 
of J = data extends out towards infinity, and so could contaminate the entire evolution. 
However, previous results comparing characteristic extraction with conventional finite radius 
extrapolation indicate that it is at most a minor effect |I8]|9|- 

Since the characteristic initial data is needed only in the far field region, linearized theory 
provides a suitable approximation. That is, it is possible to construct initial data that, at the 
linearized approximation, represents the physical situation of a gravitational field with purely 
outgoing radiation produced by sources in a central region 1 1 1 1. We first solve the case of two 
equal mass objects in circular orbit around each other in a Minkowski background as a model 
analytic problem. 

We then develop and implement a procedure to calculate characteristic metric data that 
contains purely outgoing radiation. This is done within the context of characteristic extraction, 
so that initial data on the null cone is constructed that is compatible with given data on a 
worldtube F at constant radius. We are then able to compare the waveforms computed by 
characteristic extraction using as initial data (a) J = 0, and (b) the Unearized solution. We 
find that while the dominant gravitational wave modes are largely unaffected by the choice 
of initial J, a small residual difference is visible between the two approaches, and can take 
several hundred M to be damped below other effects. 

Any mis-match between the linearized solution and the actual data is an indication of 
an ingoing radiation content. Now, on the worldtube T, the characteristic metric data is 
determined entirely by the 3+1 data so that any ingoing radiation can be traced back to an 
ingoing radiation content in the conformally flat 3+1 initial data. In this way we show that 
the 3+1 initial data contains a component of ingoing radiation which results in a long-lasting 
transient. 

The plan of the paper is as follows. Sec. |2] introduces the notation, and reviews results 
needed from previous work. Sec. [3] applies linearized characteristic theory to calculate metric 
data for two equal mass non-spinning black holes in circular orbit around each other. Sec. [4] 
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describes the method to construct metric data everywhere from data on a worldtube. Sec.|5] 
describes our numerical resuhs, which are obtained within the context of a binary black hole 
inspiral and merger. The paper ends with Sec. |6j Discussion and Conclusion. 

2. Review of results needed from other work 

2.1. The Bondi-Sachs metric 

The formalism for the numerical evolution of Einstein's equations, in null cone coordinates, 
is well known lfT2l [T3l [141 [TSl [161 [iTl . For the sake of completeness, we give a summary of 
those aspects of the formalism that will be used here. We start with coordinates based upon a 
family of outgoing null hypersurfaces. Let u label these hypersurfaces, x"^ {A — 2, 3) label 
the null rays, and r be a surface area coordinate. In the resulting — {u, r, x^) coordinates, 
the metric takes the Bondi-Sachs form |[T2l[T8l 

ds^ = - {e^P{l + Wrr) - r^liABU^U^) du" 

- 2e^f^dudr - 2r'^hABU^ dudx^ + r^hAndx^dx^, (1) 

where h^^hBc — and det{hAB) — det{qAB), with qab a metric representing a unit 2- 
sphere; Wc is a normalized variable used in the code, related to the usual Bondi-Sachs variable 
V hy V = r + WcT^. As discussed in more detail below, we represent qab by means of a 
complex dyad qA- Then, for an arbitrary Bondi-Sachs metric, Hab can be represented by its 
dyad component 

J = hABq^q''/'2, (2) 
with the spherically symmetric case characterized by J = 0. We also introduce the fields 

X=VTTjJ, U^U^qa, (3) 

as well as the (complex differential) eth operators 3 and 6 |[T9l . 

In the Bondi-Sachs framework, Einstein's equations i?^^ = 87r(rQ^ — ^g^pT) are 
classified as: hypersurface equations - Rii,q^RiA, h"^^ Rab - forming a hierarchical set 
for /?, U and Wc', evolution equation q^q^ Rab for J; and constraints Roa- An evolution 
problem is normally formulated in the region of spacetime between a timelike or null 
worldtube F and future null infinity (i7+), with (free) initial data J given on w = 0, and with 
boundary data for /3, U, U,r, Wc, J satisfying the constraints given on the inner worldtube. We 
extend the computational grid to J'^^hy compactifying the radial coordinate r by means of a 
transformation 

r 

r ^ x ~ . (4) 

r + r-wt 

In characteristic coordinates, the Einstein equations remain regular at under such a 
transformation. 

The free initial data for J essentially determines the ingoing radiation content at the 
beginning of the evolution. For the case of binary black hole evolutions in a 3h-1 formalism, 
the initial Cauchy data is determined by solving the Hamiltonian and momentum constraints, 
usually under the assumption of conformal flatness. Compatible null data initial solutions 
are not known, and so we must choose an ansatz for J which is approximately compatible. 
Previous work I8]j9| has simply set J = 0. In Section |3] we propose a refinement whereby 
J is set according to a linearized solution which is determined by the Cauchy initial data 
solution. 
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2.2. The spin-weighted formalism and the S operator 

A complex dyad is a 2-vector whose real and imaginary parts are unit vectors that are 
orthogonal to each other Further, represents the metric, and has the properties 

q^qA = 0, q^qA = 2, qAB = ^(^Ato + qAqs)- (5) 

Note that qA is not unique, up to a unitary factor: if qA represents a given 2-metric, then so 
does = e^'^qA- Thus, considerations of simplicity are used in deciding the precise form of 
dyad to represent a particular 2-metric. 

Having defined a dyad, we may construct complex quantities representing angular tensor 
components on the sphere, for example Xi = TAq^, X2 = T^^qAqs, = T^^ qAqsff ■ 
Each object has no free (angular) indices, and has associated with it a spin-weight s defined 
as the number of q factors less the number of q factors in its definition. For example, 
s{Xi) — 1, 5(^2) — 0, s{X:i) = —3, and, in general, s{X) = —s{X). We define derivative 
operators 3 and 3 acting on a quantity V with spin- weight s 

W = q^dAV + sTV, W = q'^dAV - sfV (6) 

where the spin-weights of BV and BV are s + 1 and s — 1, respectively, and where 

T = -^q^q^'WAqs. (7) 
Some commonly used dyad quantities are 



Spherical polars Stereographic 



ds2 ^ fig2 ^ sin^od(j)^ 4(c^q2 ^ ^^2)/ ^ ^2 ^ ^2^2 

q^= {l,isme) ^{1 + q^ +p^)il,i) 

T = — cot 9 q + ip. 



(8) 



(9) 



The spin-weights of the quantities used in the Bondi-Sachs metric are 

siWc) = s{/3) = 0, s(J) = 2, s( J) = -2, 

s{K)=0, s(C/) = l, s(;7) = -l. 

We will be using spin-weighted spherical harmonics ll20l I2TI gYtm, where the suffix ^ 
denotes the spin-weight, and in the case s — the s will be omitted i.e. Y^m = oYim- It 
is convenient to make use of the formalism described in ll22l fTTI . and have basis functions 
whose spin-weight components are purely real; following ifTTl . these are denoted as s^im- 
Note that the effect of the 3 operator acting on Zgra is 

3Z,„ = y£(ZTT) (lOfl) 

= V(^-lK(^+l)(^ + 2) 2^£m. 



2.3. Solutions to the linearized Einstein equations 

Ref. ifm (see also 1231 ) obtained solutions to the linearized Einstein equations in Bondi-Sachs 
form using the ansatz 

r,x^) = ^ifi^„i{r) e-Kp{ivu)) sZi^rn (H) 

for a metric coefficient F with spin- weight s. Here, we need the results for linearization 
about a Minkowski background, in which the spacetime is vacuum everywhere except on a 
spherical shell at r = tq. Strictly speaking, we should be performing the linearization about 
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a Schwarzschild (or even Kerr) background rather than about Minkowski. In the Kerr case it 
is not known how to do so, and in the Schwarzschild case the difference is that, in Eq. ( 12^? i, 
the term is replaced by a term whose leading-order behaviour is also but which is 
not representable analytically. We consider the lowest order case £ = 2 and in the exterior of 
the shell (i.e., r > tq), and describe that part of the solution that represents purely outgoing 
gavitational radiation. 

P2,v{f) = hi (constant) (12a) 

j,Ar) - (12&1 + 6^.c, + ^u-c,)^ ^ V^c, ^ 76c. ^^^^^ 

9 r 6r-^ 

= ^ + — + ^ 



18 r 
2ivci C2 \ 
~ 3r3 ~ 2^y 

W2,y[r) = r 

-66i + 12ivci + 2iv'^C2 



(12c) 



3 

2iuC2 C2 
72 



2l^^C2 



(I2d) 



r 

The solution is determined by setting the constant (real valued) parameters hi, ci and c.. The 
gravitational news corresponding to this solution is given by 

M = 3fi(n2,^ exp(wu)) 2Z2,m with n2^v = -^^^2—- (13) 

6 

We will also need the solution in the case = 0, £ = 2 in the exterior region r > tq 

P2.o{f) = bQ (constant) (14fl) 



3 

r 



2c3 




r 




2c3 


3c4 







(14^) 
(14c) 



W2M{r) = -2bor-^. (Ud) 



r 



in terms of the additional parameters 6o, C3 and C4. 



3. Black hole binaries in circular orbit: Solution in the linearized limit 

The linearized solution described in the previous section will be used to set initial data on the 
null cone. We seek a solution which corresponds roughly to the source of the gravitational 
radiation which we will eventually measure, namely a binary black hole system. In order to 
be able to apply the linearized theory, we model each black hole as having a matter density 
that is described by a Dirac-(5 function whose location moves uniformly around a spherical 
shell. More precisely, the matter density p in the spacetime is 
M 

p = —Sir ~ ro)S(9 - 7r/2) (SU - vu) + 5U -vu~ tt)) , (15) 
To 
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with respect to Bondi-Sachs coordinates (u, r, 6, (j)). The mass of each black hole is M, the 
circular orbit has radius tq, and the black holes move with angular velocity v. We next express 
p in terms of spherical harmonics 

P = ^5R(p£,mexp(|TO|wu)) (16) 

t.m 

and apply the usual procedure, that is multiplication by ^ followed by integration over 
the sphere, to determine the coefficients pt^m- We find that, for t < 2, the only nonzero 
coefficients are 

Po,o = 5{r ro)^, p,,o = -5[r - ro)|| y^, ^^^^ 

In linearized form, the Rn Einstein equation is 

P^r = 'i'Krpvl, (18) 

where vi is the covariant component of velocity in the /--direction. Imposing the gauge 
condition that the coordinates should be such that, on the worldline of the origin, the metric 
takes Minkowski form, it follows that (3 — there and consequently at all points within 
r < tq. Expanding /3 in terms of spherical harmonics 

^ = V'5R(&f,™exp(|m|wu))Zf,m, (19) 



(20) 



and integrating Eq. (18 1, we find the coefficients hi,m for r > tq. 

The determination of the remaining metric coefficients depends on the value of {£, m). 
The case (0, 0) is straightforward, and we find for r > rg 

J^O, U^O, 13^ ^, Wc = ^ + ^. (21) 

The case (2, 0) uses the results for a static shell on a Minkowski background in ifTTl . We 
solve the jump conditions across the shell for the various metric quantities ||] The result is, in 
the interior, 

b2fi = 0, j2,o(r) = , (22) 

and in the exterior 

O2.0 = , (23fl) 

The cases (2, ±2) use the results for a dynamic shell on a Minkowski background 
in ifTT I [U The script constructs the general solution inside and outside the shell r — vq, 

X Maple script for this puipose (nuO_regular_0 . map with output in nuO_regular_0 . out) are provided in the 
supplementary data. 

§ The calculation is provided in the Maple script regular_0 .map with the output in regular_0 . out in the 
supplementary data. 
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and uses the constraint equations, as well as regularity conditions at the origin and at infinity, 
to eliminate some of the unknown coefficients. It imposes the jump conditions at the shell, 
and finds a unique solution for the remaining unknowns. The result in the exterior is 

/3 = (62,2 exp(2ii/u)) ^2,2 + 3* (-«62,2 exp(2izyu)) ^2,-2 (24a) 
^ = 3fi(j2,2(r-)exp(2wM)) 2^2,2 + 3* (-ij2,2(r) exp(2z:yu)) 2^2,-2(24^) 

where 

62,2 = (25) 

ro 



and j2.2{r) takes the form given in Eq. ( I2b 1. The gravitational news is 

/fi 
6 

+ 3fi(-(2z/)3c2^ exp(i2i.u)) 2^2,-2. (26) 
6 

Although the coefficient C2 is complicated, it can be expressed to leading order in roi/ 

C2 = ^"'^o- (27) 
It follows that, again to leading order in tqI', the news is 

/27r 

(3?(-i exp(2wu)) 2^2,2 + exp(2iz/u)) 2^2,-2) (28) 

from which it is easy to deduce, via the Bondi relation, that the rate of energy loss of the 
system is 

-M^vtry'^- (29) 
o 

In the limit of a low velocity circular orbit, t^i — 1,^^ — M/ (4rQ), the above formula reduces 
to 

dE 2 

which is identical to that found from the standard quadrupole formula 



4. Constructing the metric from data on a worldtube 

In characteristic extraction, the Cauchy evolution provides the characteristic metric variables 
/3, J, U and Wc on the worldtube F, decomposed into spherical harmonics sy«,m, at every 
time step. In this section, we develop a method to find coefficients of the linearized 
solutions that provide a fit to the actual numerical data at the worldtube (to linear order and 
excluding incoming radiation). Then we use the linearized solutions with the coefficients 
just found to predict J everywhere at some chosen time u, and in this way provide initial 
data for a numerical characteristic evolution. We restrict attention to the dominant modes 
s^2,2, 8^2,-2, si^2,o- The method uses a Fourier decomposition in the time domain, and works 
well when the data is approximately sinusoidal, with amplitude and frequency varying slowly. 
Accordingly, for the binary black hole computation, the method is applied over a time domain 
that excludes both the junk radiation and the merger. 
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A metric variable A may be written as 

A = ay,2 sY2a + dYfl s5^2,0 + ay_2 sY2~2 (31) 

where * denotes the complex conjugate. The relationship between the coefficients of sl2,2 and 
s^2,-2 follows theoretically from the requirement that the spin weight metric components 
must be real; and further the metric data has been checked to confirm that it does satisfy the 
relationship. Transforming to the sZi^m basis, we find 

A — az,2 s-^2,2 + aZ,-2 s-^2,-2 + flZ.O s-^2,0, (32) 

where 

az,2 = y23?(ay,2), az -2 = -V2'^{aYa)i az,o = ay.o (33) 

so that the metric data on F can be re-expressed as coefficients of ^^2,27 s-^2,-2 and s^2,o at 
discrete time values. Although the data is oscillatory in time, it is not at constant frequency but 
is a superposition of multiple solutions with different frequencies. The linearized solutions 
behave as e"'" for fixed v, so for the theory to be applicable the next step is to decompose 
the metric data into a superposition of constant frequency components. This is achieved by 
making a discrete fast Fourier transform of each metric coefficient 

az,2,. = V exp ( (34) 



3 = 1 ^ 

where there are L data points over the time interval (ui, ml). The frequency i' is related to k 
by 

[UL - ui)L 

We found that J at from the linearized solutions provides a smoother fit to the actual 
data if high frequencies are eliminated (compare Sec.|5] Fig.[T]), and so we undertake further 
processing only for k < Li (with 1 < Li <C L); the setting of the Fourier coefficients for 
fc > ii is described later. 



In the case k > 1, for each value of k Eqs. ( 12a 1 to ( 12di evaluated at the worldtube are 
four equations for the three unknowns 61 fc, ci.fe, C2.k- Such an over-determined system can be 
tackled by a least-squares-fit algorithm, or alternatively by ignoring one of the equations so 
making the system uniquely determined. We found that the reconstructed linearized solution 



gave a better fit to the actual data at J'~^ in the case that Eq. {I2d\ for Wc was ignored. 
This then means that a comparison between the actual and reconstructed data for Wc at the 
worldtube provides an indication of the error, which is expected because of (a) incoming 
radiation in the initial Cauchy data, (b) Fourier transform effects, and (c) other effects. We 
discuss item (b) at the end of this section, and items (a) and (c) in the next section. 



In the case k = 1, u = 0, and Eqs. ( 14fli to (14^ are four equations for the constants 



60, C3, C4. We solved four equations for three unknowns using a least-squares-fit algorithm, 
because this approach led to the reconstructed data having a better fit to the actual data than 



in the case that Eq. ( 14-di was ignored. 

In this way, for a given spherical harmonic say Z2.2, we obtain values for the constants 
of the frequencies represented by fc = 1 • • • Li. Now, our purpose is to use the worldtube data 



to estimate J off the worldtube. From Eq. ( 12b 1 we can write 



j2,k(r)=do,k-\ -H 3- (36) 
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where for 2 < fc < Li 



do.k = (l2&i,fe + Givci^k + w^C2^k) , (37fl) 
= 2V6ci,fc, (37/7) 

d2,fe = ^ ci,fc, (37c) 



and for fc = 1 



do,i = = d2,i = 2^6. (38) 



We then apply the inverse discrete Fourier transform to find 



,/N 1 V- , 2m{k - l)[u - ui){L - l)\ 

Mu) = T do.k exp ^ ^ ^ ^ , (39) 

V {uL-ui)L J 

where do^k = for Li + 1 < fc < L — Li + 1, and 

do,L-k+i ^ for A: = 1 • • • (Li - 1). (40) 



Eq. (40 1 follows from the condition that cIq^u) (and all other coefficients in the time domain) 
are real. The functions c?i (u) and d2 (u) are found in a similar way, and so we are able to find 
the coefficient of J{u, r) of a given spherical harmonic, say 2^2,2- Repeating the calculation 
for the other spherical harmonics 2^2,-2 and 2^2,0 leads to a prediction of J{u,r,x^) to 
lowest order £ = 2. 

The coefficient of 2^2.0 is not oscillatory, but can rather be described as slowly varying. 
While in principle, such behaviour can be represented by a Fourier decomposition, we 
found that a better fit was obtained by regarding the solution as almost constant and solving 
Eqs. ( 14a 1 to ( 14c 1 at each time step, with Eq. ( 14-d\ used as a measure of the error[ii] 



We investigated the possibility of errors introduced in the Fourier transform and inverse 
transform process by comparing on the worldtube the actual and reconstructed values of /3, 
J and U , because the construction is such that they should be identicaj^ The following 
comparison is for the RlOO case as specified in the next section. We found that there 
was essentially no difference between the original and reconstructed data, apart from minor 
variations over about the first and last ±30 Af of the time interval (of total duration 1290M), 
presumably caused by the cut-off of high frequencies. This test was performed for both 
Li = 50 and Li = 100 with no visible difference seen in the graphs, indicating that the 
precise way in which high frequencies are removed is not important. 



5. Numerical results 



We apply the method outlined in the previous section to the problem of measuring 
gravitational waves from binary black hole simulations. Following our implementation of 
characteristic extraction |i8n9J, we first evolve a spacetime with a 3+1 (Cauchy) evolution 
code, recording metric data on a world tube, F, of fixed coordinate radius. This data is 
subsequently used as inner boundary data for a null-cone evolution of the Einstein equations, 
which transports the data to J'^, where the gravitational waves are measured. The linearized 

II The Matlab scripts f t_wt_driver . m, FT_WT . m and nuO . m are provided in the supplementary data. 

^ Note that Wc is not necessarily identical, since there are may be differences due to nonlineaiities and incoming 

radiation. 
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Worldtube location Initial time Initial data 



Data set 



R.r{M] ua[M] J 



JO-RlOO-uO 
J0-R250-U0 



100 
250 
100 
100 
100 
250 






450 
900 
450 
900 



J ^0 
J = 
J = 
J = 
J — Jli, 
J = Jli, 



J0-R100-U450 
J0-R250-U900 
Jlin-R100-u450 
Jlin-R250-u900 



lin 



lin 



Table 1. The characteristic evolutions performed based on the same Cauchy evolution of two 
equal mass non-spinning black holes. 



solution allows us to specify data for J on the initial null cone which is compatible (to the 
linear level) with the data on F and thus, importantly, to the Cauchy 3+1 spacetime. 

As a fiducial test case, we return to the well-studied model of an 8-orbit binary system 
with equal mass non-spinning black holes carried out in |[8]|9l. For the Cauchy evolution, 
we use the Llama multipatch code described in ||25] |251 . We output metric data on two 
worldtubes located at i?r = lOOA/ and i?r = 250M that are used as inner boundary data 
for a subsequent characteristic evolution. The waveforms at should be independent of the 
worldtube location. Thus, evolutions from different worldtube locations help us validate our 
results. Table [T] summarizes the various characteristic evolutions that we have performed, all 
of which are based on the same Cauchy data, but with different characteristic initial data, J, 
and starting points in Bondi time, uq. 

The first approach follows the original prescription laid out in |[8l|9l. The characteristic 
evolution is started coincident with the first available Cauchy data, at coordinate time tg = 
uo = 0. The characteristic variable J is initialized by the shear-free solution J = 0. Since we 
are beginning the characteristic evolution from the initial Cauchy slice, these models include 
the spurious junk radiation contained in the conformally flat constraint solution. Models 
JO-RlOO-uO and J0-R250-u0 listed in Table [T] follow this prescription, using data from the 
worldtubes at i?r = lOOAf and i?r = 250M, respectively. 

It is interesting to compare the results of the fully nonlinear characteristic Einstein 
evolution with the corresponding linearized solution. Fig. [T| plots the {£,m) = (2,2) 
spherical harmonic modes of Jnum. computed by the null Einstein evolution JO-RlOO-uO. The 
linearized solution Jiin is computed using linearly reconstructed worldtube data according to 
the prescription in Sections [3] and |4] We compute the linearized solution from the boundary 
data at F after the initial data junk radiation has passed the worldtube radius i?r, at around 
u = 150 A/, when the system has settled to the expected binary black hole inspiral pattern 
compatible with the solution construction. The upper panel of Fig. [T] plots the real and 
imaginary parts of J, evaluated at The center panel plots the amplitudes of J, while the 
bottom panel shows the relative difference between the linearly estimated and the fully 
relativistic result J„ui„ (model JO-RlOO-uO). The linearized Jun and numerically evolved 
Jnum differ initially, but eventually, after around u — 450 A/, differ by less that 1%, which 
remains consistent for the bulk of the time. 

We first discuss whether the difference between the solutions can be due to effects other 
than ingoing radiation; such effects could include (a) nonlinearity, or (b) the linearization 
background being Minkowski rather than Schwarzschild. Looking at Figs.[T]and|2] we see that 
to lowest order the metric components are slowly varying sinusoidal functions; this statement 
also applies to the other metric components (graphs not shown). We would therefore expect 




Figure 1. Components of the 2^^2,2 mode of J at J''^ estimated from boundary data at 
_Rr = lOOA/ once using linearized solutions, and once using the full non-linear characteristic 
evolution JO-RlOO-uO. The numerical solution is denoted by Jnum, while the reconstructed 
hnear solution is J^^ . The linearized solution makes use of data starting from a time after 
which the initial burst of junk radiation has left the system. The two solutions agree reasonably 
well only after a time i2 Eq. j4TJ, a time after which the incoming radiation content of the 
Cauchy initial data has essentially settled to zero. 



that the magnitude of effects (a) and (b) would be roughly constant, if not with some increase 
at later times as merger is approached. Indeed for nonlinear effects we can look at the Einstein 
equation for Rn, which is /3 = to linear order, with the actual value being an indication 
of the magnitude of nonlinearity, and we find that this quantity does slowly increase with 
time. However, Figs. [T] and |2] show that the difference between the solutions, until about 
450M, is getting smaller. This indicates that the effect is independent of non-linearities 
and present already at the linear level. The linearized solution assumes that the radiation is 
purely outgoing, whereas the actual data may contain incoming modes originating in either the 
characteristic or Cauchy initial data - both options being possible since J at is influenced 
by both data sets. Thus, the explanation for Fig.[T]is that it reflects the slow decay of the effect 
of incoming modes in the initial data, until saturated by other factors (such as non-linear 
effects). 

A similar effect is seen in the characteristic variable Wc, related to the Newtonian 
potential, plotted in Fig. |2] In this case, however, there is clarity about the source of the 
incoming radiation: it must be in the Cauchy data. This is because in characteristic extraction. 
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Figure 2. Components of the 2^2,2 mode of Wc at the worldtube ijp = lOOAf once using 
the linearly reconstructed data, and once using the original data as obtained from the Cauchy 
data. The numerically obtained Cauchy data is denoted by Wc, num. while the reconstructed 
data is VF^, nn. The reconstructed data makes use of the original Cauchy data starting from 
a time after which the initial burst of junk radiation has left the system. The numerical and 
reconstructed data agree reasonably well only after a time t2 Eq. {41) , a time after which the 
incoming radiation content of the Cauchy initial data has essentially settled to zero. 



the characteristic metric at the worldtube is determined entirely by the Cauchy data. Again, 
the lower panel shows an approximately exponential decay in the differences, until around 
u = 400Af . The residual steady state differences result from other effects, which gradually 
increase with the strength of the gravitational radiation towards the binary merger. 

The findings above indicate that a physically expected purely outgoing inspiral radiation 
pattern is only present after some time 

U > ^incoming ; (41) 

where ^incoming is the length of a time interval until the incoming radiation content of 
the Cauchy initial data has settled to a negligible amount at the given worldtube location. 
Hence, in order to construct physically meaningful and consistent initial data via the outgoing 
linearized solution, it is optimal to begin the solution at a time uq > ^incoming after which 
both the junk and incoming radiation content of the initial data have subsided. The results of 
Figs.[T]and|2]suggest that the linearized solution provides a reasonably good approximation to 
the data u « 450Af , at which time the system has settled into an outgoing radiative solution. 
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For instance, in model Jlin-R100-u450, the exact worldtube location is Rr = 100.8492. 
Allowing 50Af for the visible outgoing junk radiation to pass, we set the time range 
over which we use the boundary data for linearized solution construction to {uQ,Uf ) — 
(150.192, 1439.856), which includes the inspiral, but not the merger and ringdown. The time 



increment in is given by du = 0.144, thus comprising 8967 data points. Referring to Eq. ( 34 1, 



we have L = 8957, Li = 100. In order to determine J^n at the initial time uq, we use the 



general form of the linearized solution, Eq. (2.3 i 



^ - (eo + - + ^) 2Y2,2 + (eo + - + ^) 2^^2,-2 
\ r / \ r / 

+ fea + - + ^) 2Y2.0 (42) 

The coefficients are determined by comparing with the worldtube data. We perform a Fourier 
transform on the numerically determined worldtube variables over the time interval (uq, u/). 



The spectrum determines the constants 61, ci and C2 of Eqs. 2.3 at each fixed frequency v. 



These values are transformed back to the time domain, and evaluated at t = 450A/ in order 



to determine the coefficients of Eq. (42 1 



eo = (4.5217 + 3.7702i) x 10-^ ei -0.04578 - 0.17159i, 

62 = 12.582 + 42.500i, 63 = -1.4788 x lO"^, (43) 

64 = 0.020365, 65 = -35.563. 

A goal of this paper is to investigate, within the context of characteristic extraction, 
the effect of the initial data on the calculation of the gravitational news. To this end, we 
compare waveforms at from two characteristic evolutions based on the same Cauchy 
boundary data, but different initial data constructions: Jlin-R100-u450 and J0-R100-u450. 
Both evolutions use boundary data from Ry = lOOM and begin at the initial time uq — 
t2 = 450M, which was determined above to be a point where the linearized solution is 
well-matched to the nonlinear evolution. The model Jlin-R100-u450 uses the initial data 
determined by the linearized solutions, Eqs. (42 43}. In contrast, J0-R100-u450 simply sets 



J = 0, corresponding to the original prescription of ||8]|9l. Fig.plplots the Bondi news at J'^ 
as computed from both evolutions, denoted by A/{f^° and Ap'^'^or the linearized and J = 
initial data runs, respectively. Whereas the phase of Af shows very little difference between 
the runs (middle panel), the amplitude shows visible oscillations for the A/jf evolution (upper 
panel and inset). The waveforms agree to within 1% only after a time u = 400M (which must 
be added to the uq — 450A/ starting point of the simulation). Therefore, the influence of the 
J = ansatz for the characteristic initial data has a notable impact over an extended time. 

We note, however, that the original prescription of 1,8, |9J used characteristic initial data 
J = at the initial Cauchy time uq = to — Q (that is, including the junk radiation). At that 
time, the shear-free approximation J = 0, for the characteristic initial data is compatible with 
the Cauchy initial data solution. 

Hence, we also compare the waveforms of model Jlin-R100-u450 against those of model 
JO-RlOO-uO, where the latter model uses J = initial data at time uq — 0. The results are 
plotted in Fig. |4] with the JO-RlOO-uO results labeled by A/J. We still observe an oscillation 
in the amplitude in Ao, though it is drastically reduced compared to the A/'q of model JO- 
R100-u450 shown in Fig. [3] The relative errors in amplitude, plotted in the bottom panel of 
Fig. [4] are well below 1% over the entire evolution, and the total dephasing is smaller than 
A(j> = 0.04rad. We note that the differences are larger than the systematic error reported in 
Is, 9|. In that work, the error estimate refers to the difference between evolutions starting 
from the same initial data but different worldtube locations, at a fixed resolution (though the 
results converge to the same waveform as the resolution is increased). Our results indicate that 




Figure 3. Time domain differences between the Mq^^ and waveforms of models JO- 

R100-U450 and Jlin-R100-u450, computed from worldtube location i?r = lOOA/ for which 
the characteristic runs have been initialized by J = and linearized solutions at a time 
«0 = 450J\/, respectively. The top panel plots the wave amplitude, the middle shows the 
phase, and the lower plots the differences between the two solutions. The A/^f 50 data shows 
notable oscillations in amplitude at early times (inset), which decay exponentially with time. 
The waveforms have been aligned at the amplitude peak. 



at typical current resolutions, the choice of characteristic initial data has a larger influence on 
the simulation error than the worldtube location. The initial burst of junk as well as spurious 
incoming transients can alter the measured waveforms by an amount that is of the order of 
the discretization error of current numerical relativity codes (e.g. [21 \ and references therein), 
over a period of several hundred M. 

The numerical tests described so far were with the extraction radius i?r = lOOA/. All 
these runs were repeated with the extraction radius re-set to i?r = 250A/. The results are 
qualitatively similar to the i?r = lOOM case, and the details are not presented here. One 
interesting feature was the behaviour of Wc at the worldtube, i.e., the analogy to Fig. |2] it 
took until about 900Af until the decay in the difference between the linearized and actual 
data was saturated. This is surprising since one would expect that the radiation that passed 
i?r = lOOA/ at about u = 450A/ would pass i?r = 250A/ at w = 600A/. One possible 
explanation is that the saturation is due to nonlinear effects, and since they are somewhat 
weaker at i?r ~ 250Af, saturation takes longer. Furthermore, the boundary data amplitudes 
are one order of magnitude smaller at i?r = 250A/ than those at Ry = lOOAf , and thus 
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Figure 4. Time domain differences between the A/q (starting off at mq = <o = 0) and ■N'y^ 
(starting off at uq = 450) waveforms of models JO-RlOO-uO and Jlin-R100-u450, computed 
from worldtube location _Rp = 100-A/ for which the characteristic runs have been initialized 
by J = and linearized solutions respectively. The top panel plots the wave amplitude, the 
middle shows the phase, and the lower plots the differences between the two solutions. The 
data shows slight oscillations in amplitude at early times (inset), which decay exponentially 
with time. The waveforms have been aligned at the amplitude peak. 

potentially more sensitive to incoming modes. Further studies would be required to fully 
understand the nature of these effects. 

6. Discussion and Conclusion 

Characteristic evolutions provide a means of determining radiated gravitational energy which 
is free from the ambiguities associated with local measures, namely non-linear effects in the 
near-zone, as well as ambiguities due to gauge and extrapolations to infinity. The principal 
remaining issue has been to specify appropriate initial data on the null cone, compatible 
with the data on the worldtube and the Cauchy 3-1-1 spacetime. The strong correlation 
between finite radius results and characteristic extraction, as well the invariance of the results 
on the worldtube location observed in IS] |9] suggests that even the simplest initial data 
ansatz, J ~ 0, can provide results which are accurate enough for astrophysical estimates, 
provided the Cauchy 3-1-1 spacetime is essentially free of radiation at the initial characteristic 
time. However, in scenarios where strong outgoing radiation is present during the initial 
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characteristic time, J = data effectively represents incoming radiation, and may alter 
significantly the evolution of the wave signal towards 

The gravitational wave solution developed in Sec. [3]provides initial data J — which 
is compatible, to the linearized level, with outgoing radiation from a binary system. As such, 
it provides a more physically motivated starting point than the shear-free, J = 0, alternative. 
Importantly, we find that the evolutions which take place from either J — Jun initialized 
at a time uq = 450Af when outgoing radiation is present at the worldtube, and J — 
at the initial time Ug = io = when the Cauchy 3+1 spacetime is conformally flat are 
very similar (compare Fig. |4|l, and indeed very similar to the purely 3+1 result which can be 
obtained by polynomial extrapolating finite radius measurements. That is, for simple choices 
of initial J, the physical conclusions are not altered dramatically, provided the Cauchy 3+1 
spacetime does not contain strong amounts of radiation at the worldtube location during the 
initial characteristic time uq. 

On the other hand, we have demonstrated that the choice of characteristic initial data 
does result in a small but measurable difference, which decays at a slow exponential rate 
over a time period of several hundred M. We see this transient in the graphs of gravitational 
news exhibited in Fig. [3] Since the linearized initial data contains only an outgoing mode, we 
conclude that the shear-free characteristic initial data contains incoming radiation. While this 
is expected, it is interesting that it takes so much time for the effect to decay away. 

This study was designed to assess the long-term effect of characteristic initial data, but 
it also provides information about the incoming radiation in 3+1 initial data. The point is that 
the characteristic data on the worldtube is determined entirely by the 3+1 data, and the extent 
to which this data does not fit the linearized solution is a measure of its incoming radiation 
content. By construction, the quantities /?, U and J in the linearized solution must fit the data, 
with the difference in Wc being an indication of incoming radiation in the 3+1 evolution. 
Since Wc is not a gauge invariant quantity, it is not possible to make a quantitative statement 
about the magnitude of the incoming radiation, but Fig.[2]indicates that it takes until at least 
±400Af until it is possible to neglect the effect of incoming radiation in the region r < lOOM. 

The effect of incoming radiation in both 3+1 and characteristic initial data decays 
exponentially, but even so it has a more long-term impact than the outgoing junk radiation 
which passes the r = lOQM worldtube radius by approximately u = 15QM. The effect 
of incoming radiation has not received much attention in previous literature dealing with the 
problem of initial data construction. It has potential significance for the construction of high- 
accuracy gravitational-wave templates. Further investigations may also explain properties of 
the incoming radiation content, giving further insight to the peeling property of more complex 
spacetimes other than Kerr. In particular, it will be important to determine, using a gauge 
invariant measure, the long-term effect of incoming radiation in 3+1 initial data sets. 

The results suggest some promising avenues for future study. Current methods 
("characteristic extraction") transport data in one direction (from the Cauchy to the 
characteristic code). Great efficiencies would be possible if the coupling were also earned out 
in the other direction, so that the characteristic evolution would provide outer boundary data 
for the Cauchy evolution. The linearized wave solution provides a simple recipe for isolating 
ingoing vs. outgoing modes on the characteristic grid, and may be useful in designing a 
method for stably transporting data in both directions across the world tube interface. 
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